Beyond Eliashberg superconductivity in MgB2: anharmonicity, two-phonon scattering, 

and multiple gaps 
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Density-functional calculations of the phonon spectrum and electron-phonon coupling in MgB2 are 
presented. The Egg phonons, which involve in-plane B displacements, couple strongly to the px,y 
electronic bands. The isotropic electron-phonon coupling constant is calculated to be about 0.8. 
Allowing for different order parameters in different bands, the superconducting A in the clean limit 
is calculated to be significantly larger. The E2g phonons are strongly anharmonic, and the non-linear 
contribution to the coupling between the -B29 modes and the px,y bands is significant. 
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The recent discovery of superconductivity near 40 K 
in MgB2 has generated much interest in the properties of 
this simple intermetallic compound Q. A significant B 
isotope effect strongly suggests phonon-mediated pairing 
1^. To explain the large Tc, an electron-phonon coupling 
(EPC) constant of A « 1 is needed. Yet estimates of 
the coupling strength based on the latest measurements 
of the low- temperature specific heat combined with 
the density of states (DOS) from density-functional cal- 
culations [Q, yield A « 0.6-0.7. Further, the measured 
temperature dependence of the electrical resistivity ^ is 
consistent with Xtr ^ 0.6. First-principles calculations 
of the EPC give A « 0.7 - 0.9 Clearly there is a 

problem in reconciling all these numbers. Another puz- 
zle involves tunneling measurements of the gap. Values of 
2A/kBTc ranging from 1.2 to 4 have been reported. The 
values below the BCS weak coupling limit of 3.5 have 
been attributed to surface effects, but the best-quality 
spectra [|j show a very clean gap with 2A/kBTc = 1.25. 
Sharvin contact measurements Q reveal a gap at 4.3 
meV (2A/A:bTc = 2.6), and additional structures at 
iA/ksTc — 1.5 and 3, raising the possibility of multiple 
gaps. Careful analysis of the temperature and magnetic 
field dependence of the specific heat suggests anisotropic 
or multiple gap structure as well Thus, even if su- 
perconductivity in MgB2 is phonon-mediated, it is likely 
that an analysis beyond the simple isotropic Eliashberg 
model is needed. 

The MgB2 lattice consists of two parallel systems of 
flat layers. One layer contains B atoms in a honeycomb 
lattice, the other Mg atoms in a triangular lattice halfway 
between the B layers. First-principles calculations [Q find 
that the electronic states near the Fermi level are primar- 
ily B in character and the Fermi surface (FS) comprises 
four sheets: two nearly cylindrical hole sheets about the 
F-A line arising from quasi-2D p^.y B bands, and two 
tubular networks arising from 3D pz bonding and anti- 
bonding bands B. The difference in character between 



the sheets raises the possibility that each has a distinct 
gap that could be observed in the clean limit. Such in- 
terband anisotropy enhances the effective EPC constant 
relevant to superconductivity and decreases the coupling 
constant for transport, compared to the average values 
pO|-p^. This could explain the discrepant values of A 
deduced from different types of experiments. 

In this paper, we report first-principles calculations of 
the EPC in MgB2. The coupling constant is decomposed 
into contributions from the four different bands crossing 
the Fermi level, allowing for an analysis of the effects 
of interband anisotropy on the superconducting Tc and 
gap structure. The strongest coupling arises from the 
E2g phonon modes along the F to A line, which strongly 
interact with the quasi-2D electronic states. This phonon 
mode is calculated to be highly anharmonic, and it also 
has significant nonlinear contributions to the EPC psf . 

Harmonic phonon frequencies and linear EPC param- 
eters were calculated using the linear-response method 
within the local density approximation Jl^ . Norm- 
conserving pseudopotentials were used, with a plane- 
wave cutoff of 50 Ry. We used the experimental crystal 
structure, with a = 3.08 A and c/a~ 1.14 The elec- 
tronic states were sampled on grids of up to 24'^ k-points 
in the full Brillouin zone, and the dynamical matrix was 
calculated on a grid of %^ phonon wavevectors q p^ . 

The calculated phonon density of states F{uj) and 
Eliashberg function a^F(w) are plotted in Fig. |^. The 
results are similar to those reported in Refs. Q and [Q. 
All of these calculations give a slightly softer phonon 
spectrum than what is observed in neutron experiments 
llTsl . While F and a^F are similar in shape in many 
materials, they are strikingly different in MgB2. In 
particular, a^F has a pronounced peak in the range 
of 60 to 70 meV arising from dispersive optic modes 
that do not give rise to large structures in F. Corre- 
spondingly, the average phonon frequency ujave — 55.3 
meV is less than the logarithmically averaged frequency 
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Win = exp[A~-^ J In uja'^F{cu)uj~^du}] = 56.2 meV, despite 
the fact that logarithmic averaging preferentiahy weights 
lower frequencies. The isotropic EPC constant, which 
determines Tc in the dirty limit, A^^ = 2 / LO~^a'^F{u})dcj 
is found to be 0.77, in reasonable agreement with other 
calculations ^,^,0- 

The peak in a F between 60 and 70 meV arises from 
the E2g phonon modes with q along the F-A line. This 
Raman active phonon mode, doubly degenerate at F, in- 
volves in-plane, hexagon-distorting displacements of the 
B atoms. In fact, by symmetry, this is the only mode 
at F that has a linear EPC. Going away from the F- 
A line the EPC drops sharply when the phonon wave 
vector q becomes larger than the diameter of the 2D 
Fermi surface; at the same time the frequency increases 
by roughly 30%. This indicates that the reason why this 
B-B bond-stretching mode is not the highest-frequency 
mode at F is because of softening due to EPC. How- 
ever, this softening should weaken in the superconducting 
state, since some of the screening electrons form Cooper 
pairs and are removed from the Fermi sea. Such hard- 
ening of optical phonons below Tc was first discussed by 
Zeyher and Zwicknagl [|l9| . The overall scale of the rela- 
tive hardening, Auj/lo, is set by a specific EPC constant, 
Xzz = 2w"^ |fffc,fcP^(eki),where g is the EPC matrix 
element. (The Fermi level is set to zero.) In the BCS 
limit, Auj/uj is a known analytical function of uj. We 
calculate Xzz — 0.6, for the £'23 mode. Taking A ~ 5 
meV we predict about a 12% hardening of this mode be- 
low Tc- This shift should be easily observable in Raman 
or neutron experiments. 

Since the 2D FSs are calculated to play an impor- 
tant role in the EPC, we have decomposed the rele- 
vant electronic characteristics in terms of the four sheets 
of the FS. We Hst in Table I the partial DOS iV, = 
and plasma frequencies ^^^^ = -^^Wi = 

I]k"k,:.Q^(<^ki)7 with i = 1(2) referring to the 
light (heavy)-hole 2D sheets of the Fermi surface, and 
i = 3(4) to the pz bonding (antibonding) sheets. The 
EPC constant was also decomposed into contributions 
from scattering of an electron from band i to band j: 

ij i 

U.,N,N, = 2E^qJlCk+q,f'^(ekO'5(ek+qj)- 
kqiy 

Here ujqi, is the frequency of the corresponding phonon, 
and is the standard (Eliashberg) isotropic coupling 
constant. Allowing for interband anisotropy of the or- 
der parameter (clean limit), the effective coupling con- 
stant for superconductivity A^^^ is given by the max- 
imum eigenvalue of the matrix Ay = UijNi, which is 
always larger than A^^. Assuming the same interac- 
tion parameters Uij for transport properties, the low- 
est order variational approximation for the Boltzmann 



equation corresponds to the transport EPC constant 
^tr = Yli^i^i/^ ■ Oil the other hand, allowing vari- 
ational freedom for the different sheets of the Fermi sur- 
face yields an effective transport coupling constant which 
is always smaller than A^^. In effect, the different bands 
provide parallel channels for conduction, so that when 
"scattering-in" is neglected, W/X^/J = Y.^W^/X^ 

The calculated interaction parameters Uij are listed in 
Table II. Due to similarities between the two 2D sheets, 
and between the two 3D sheets, we have simplified the 
model to allow for two different order parameters for 
these two sets of bands. This gives Uaa — 0.47 Ry, 
Ubb ^ 0.10 Ry, and Uab = 0.08 Ry, where A and 
B stand for the 2D and 3D bands, respectively. Then 
Xa = 1.19 and Xb = 0.45, suggesting de Haas-van 
Alphen mass renormalizations of ~ 2.2 and ~ 1.5, for 
the two sets of bands, and specific heat renormalization 
of 1.77 The resulting anisotropic effective coupHng 
constant for superconductivity is A^^-'' = 1.01. Using 
the Allen-Dynes approximate formula for Tc |^ , we find 
that to have Tc = 40 K, a Coulomb pseudopotential of 
/i* « 0.13 is needed. This is more physically plausible 
than the ^* « 0.04 required when A^c is used. For trans- 
port, interband anisotropy reduces the in-plane coupling 
constant Xx,y from 0.70 to 0.57, but has essentially no 
effect on the out-of-plane A^ = 0.46 (Table HI). This 
is because the anisotropic formula accounts for the fact 
that the transport is mostly due to the 3D bands in any 
direction, simply because they couple less with phonons. 
The measured resistivity Q can be fit remarkably well 
with the Bloch-Griineisen formula using the calculated 
isotropic o^yF , with the in-plane and out-of-plane contri- 
butions appropriately averaged for polycrystalline sam- 
ples p3[ |. However, with band anisotropy, the resistivity 
is slightly underestimated. 

The temperature dependence of the individual gaps 
Ai in the weak-coupling multigap model is defined by 

A,; 



E,, C/.,iV,A, / di? tanh(^^5^)/yF" 



A2. As 

shown in Fig. ||, the larger 2D gap is calculated to be 
BCS-like, with a sHghtly enhanced 2A/Tc, while the 3D 
gap is about three times smaller in magnitude. Thus, 
in the clean limit, MgB2 should have two very different 
order parameters, which in turn should affect thermo- 
dynamic properties in the superconducting state. Ex- 
periments indicate that the coherence length in MgB2 
is close to 50 A. The mean free path corresponding to 
the residual resistivity observed in Ref. |^ is more than 
1000 A, so that 2it(^/1 « 1/3. This is in the reasonably 
clean regime, and it is likely that the intrinsic resistiv- 
ity is even smaller. However, stronger defect scattering 
should be detrimental to superconductivity: using the 
Allen-Dynes formula with the same /i* — 0.13, we get an 
isotropic Tc — 22 K. Indeed, irradiation has been found 
to drastically reduce Tc |24|. Some of the experimental 



manifestations of multigap superconductivity would be 
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a reduced and impurity-sensitive specific-heat jump at 
Tc, a deviation of the critical-field temperature depen- 
dence from the Hohenberg-Werthamer formula, a reduc- 
tion of the Hebel-Slichter peak in NMR, and a substantial 
difference between the in-plane and out-of-plane tunnel- 
ing spectra. In particular, the latter should see only the 
smaller gap | p5[ |. 

Note that A ~ 1 is in the intermediate coupling regime. 
Furthermore, the multigap scenario suggests particular 
sensitivity to impurity scattering. This means one should 
really solve the anisotropic Eliashberg equations with im- 
purity scattering, rather than the weak coupling BCS 
equations we used. Thus we do not make any quantita- 
tive thermodynamic and spectroscopic predictions here. 

We focus now on the E2g phonon modes, which carry 
the lion share of the EPC. We have examined this mode 
at r in detail using frozen-phonon calculations. The 
calculations were carried out using a general potential 
LAPW code with the same setup as in Ref. This 
mode has substantial anharmonicity. A fit of the to- 
tal energy for B displacements u between ±0.07 a.u. to 
a fourth-order polynomial {Etot — '^o.nU^) gives a2 = 
0.434 Ry/a.u.2, 03 = -1.868 Ry/a.u.^ and = 14.926 
Ry/a.u.*. Anharmonicity increases the frequency of this 
mode by about 15%, which should result in an overall 
reduction of A by ^--^10%, and an increase of win by ^6%. 

More interestingly, the i?2g modes have a signifi- 
cant nonlinear coupling with electrons. The linear cou- 
pling vertex, gi, corresponding to scattering by a single 
phonon, is proportional to matrix elements of dV/dQ, 
where Q = v2Mcju, while the second-order coupling, in- 
volving exchange of two phonons, is proportional to ma- 
trix elements of d^V/dQ^. At F, the Hellman-Feynman 
theorem allows the calculation of gi via deformation 
potentials. This is no longer the case for (72- One 
can use d^e^/dQ^ only as a qualitative estimate of 
{d^V / dQ"^) . For the cylindrical sheets of the Fermi sur- 

1/2 

face, ((c?^ek/c?Q^)^) — 26 and 20 mRy as compared 

1/2 

to {{dck/ dQ)'^) — 14 and 16 mRy. This suggests that 
nonlinear pairing via two-phonon exchange is compara- 
ble to or even larger than the linear coupling [ p6[ . The 
reason for this anomalous behavior lies in the specifics of 
the band structure of the 2D px,y bands. In the nearest- 
neighbor tight-binding approximation it can be described 
as 

4u^ = {tl+tl){6 + J2 cos G, ) + 6UU XI cos 

i i 

(Avk)^ = 4(4 - tlfiY. '^o^' cos G, cos Gj) 

i i^j 

+ 3(t„-t,)4(;XsinG,)', 

i 

where Gi = a^k, and are the three smallest lattice 



vectors. At the F point Uk = 0, and there are two doubly 
degenerate states. The bonding pair forms the 2D Fermi 
surfaces, and it appears that the main effect of the 
phonons is to lift the degeneracy at F, thereby changing 
the overall splitting between the two subbands. This ef- 
fect does not depend on the sign of the ionic displacement 
(the degeneracy is lifted either way) and thus is nonlin- 
ear by definition. This is illustrated in Fig. |3| by the two 
Fermi surface plots for two opposite E2g phonon patterns, 
both corresponding to Q « 1. One can see that while for 
the 3D sheets the coupling is mostly linear (changes of 
the Fermi surface are opposite), for the 2D cylinders it 
is mostly quadratic (changes are the same, cf. the undis- 
torted Fermi surface in Ref. Q|). Nonlinear EPC is also 
a likely source of anharmonicity: as discussed above, the 
contribution to the E2g phonon self-energy from the EPC 
with the 2D FSs amounts to as much as 20 mcV, as evi- 
denced by the softening of this phonon at F. A sizeable 
part of this softening probably comes from two-phonon 
processes. Quartic anharmonicity of the ion-ion interac- 
tion arises in the fourth order in the linear interaction 
constant, but in the second order in the nonlinear one. 

In summary, we have presented a first-principles inves- 
tigation of the electron-phonon coupling in MgB2. Inter- 
band anisotropy enhances the coupling constant from its 
isotropic dirty-hmit value of A^^ = 0.77 to an effective 
clean-limit value of X^l^ — 1.01 for superconductivity. 
With Win — 56.2 meV, this A^^-'^ is arguably consistent 
with the measured Tc of nearly 40 K. In the clean limit, 
we predict two different superconducting order parame- 
ters: a larger one on the 2D FSs and a smaller one (by 
approximately one third) on the 3D FSs. Since current 
experiments suggest that MgB2 is indeed in the clean 
limit, multiple gaps should be observable. There are hints 
of this in both the tunneling and thermodynamic data. 

The E2g phonon mode involving in-plane B motion 
provides the strongest coupling. We predict a harden- 
ing of ^ 12% of this mode at the zone center below Tc- 
In addition, this mode is highly anharmonic and it may 
also have significant nonlinear electron-phonon coupling. 
The former likely reduces the linear EPC, while the lat- 
ter increases the total EPC. Further work is needed to 
better elucidate the effect of the anharmonicity and non- 
linear coupling on the superconducting properties of this 
material. 
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FIG. 1. Phonon density of states and Eliashberg function. 




FIG. 2. Superconducting order parameters Ai in the multi- 
gap weak-coupling approximation (solid lines), and in the 
isotropic (dirty) BCS limit (thin line). The BCS order pa- 
rameter corresponding to the the same Tc as the multigap 
model is shown by the dashed line. 




FIG. 3. Fermi surfaces for two frozen-phonon patterns of 
E2g symmetry with B displacements of ±0.07 a.u. One of the 
cylindrical sheets has shrunk and is barely visible. 

TABLE I. Band-decomposition of the electronic density 
of states at the Fermi level and in-plane and out-of-plane 
plasma frequencies. The density of states is in units of 
states/Ry/spin/cell, and the plasma frequency is in eV. 





total 
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N{Ef) 


4.83 


0.66 


1.38 


1.26 1.52 




7.21 


2.91 


2.95 


3.05 5.04 




6.87 


0.44 


0.52 


4.62 5.06 


TABLE II 


Band-decomposition of the electron-phonon lu- 


teraction. 












ll 


12 


13 


14 22 


(Ry) 


0.676 


0.419 


0.064 


0.096 0.477 


ij 


23 


24 


33 


34 44 


Uij (Ry) 


0.064 


0.097 


0.113 


0.106 0.092 



TABLE III. Calculated superconducting and transport 
electron-phonon coupling parameters in both the isotropic 
limit and with interband anisotropy. The last column con- 
tains an average Xtr appropriate for polycrystalline samples. 











^tr,ave 


isotropic 


0.77 


0.70 


0.46 


0.58 


multi-gap 


1.01 


0.58 


0.46 


0.53 
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